# To Do 1. Work with a stats student (Yayie Duan) to write a Shiny app.
require(data.table)Loading required package: data.table
Registered S3 method overwritten by 'data.table':
method from
print.data.table
data.table 1.13.6 using 2 threads (see ?getDTthreads). Latest news: r-datatable.com
require(ggplot2)Loading required package: ggplot2
require(plotly)Loading required package: plotly
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching package: 㤼㸱plotly㤼㸲
The following object is masked from 㤼㸱package:ggplot2㤼㸲:
last_plot
The following object is masked from 㤼㸱package:stats㤼㸲:
filter
The following object is masked from 㤼㸱package:graphics㤼㸲:
layout
require(deming)Loading required package: deming
require(MethComp)Loading required package: MethComp
ATTENTION! In ggplot, must explicitly set ‘inherit.aes = FALSE’, otherwise ‘geom_ribbon’ expects values for aesthetics from previous geoms (i.e. tt variables ps, pb, etc.). See solution here:
(geom_ribbon expects an unused y-axis variable to be present)[https://github.com/tidyverse/ggplot2/issues/3364]
load("data/tt.RData")# Log2-transformed variables
tt[, c("l2pb",
"l2ps") := list(log2(pb), #xx
log2(ps))] #yy
# Regression on log2-transformed variables
m2 <- lm(l2ps ~ l2pb,
data = tt)
summary(m2)
Call:
lm(formula = l2ps ~ l2pb, data = tt)
Residuals:
Min 1Q Median 3Q Max
-0.67656 -0.50577 0.02698 0.40965 0.72161
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.5353 1.4238 -3.185 0.00973 **
l2pb 1.4596 0.1978 7.379 2.37e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.553 on 10 degrees of freedom
Multiple R-squared: 0.8449, Adjusted R-squared: 0.8293
F-statistic: 54.45 on 1 and 10 DF, p-value: 2.371e-05
new_l2pb <- seq(from = 5.5,
to = 8.75,
length.out = nrow(tt))
# new_l2pb <- seq(from = 5,
# to = 9,
# length.out = nrow(tt))
pi_20_80 <- data.table(predict(m2,
newdata = list(l2pb = new_l2pb),
interval = "prediction",
level = 0.6))
pi_20_80$l2pb <- new_l2pb
head(pi_20_80)p1 <- ggplot(tt,
aes(x = l2pb,
y = l2ps,
fill = lab)) +
geom_errorbar(aes(ymin = log2(10^4*(prob.stroke - 1.96*se.stroke)),
ymax = log2(10^4*(prob.stroke + 1.96*se.stroke))),
size = 0.5) +
geom_errorbar(aes(xmin = log2(10^4*(prob.bleed - 1.96*se.bleed)),
xmax = log2(10^4*(prob.bleed + 1.96*se.bleed))),
size = 0.5) +
geom_point(shape = 21,
aes(size = sqrt(N))) +
geom_abline(slope = m2$coefficients[2],
intercept = m2$coefficients[1],
linetype = "dashed") +
geom_ribbon(data = pi_20_80,
aes(x = l2pb,
ymin = lwr,
ymax = 9.25),
fill = "green",
alpha = 0.2,
inherit.aes = FALSE) +
scale_x_continuous("Log2 of Bleeding per 10^4 patients",
# limits = c(5.5, 8.5),
expand = c(0, 0)) +
scale_y_continuous("Log2 of Stroke per 10^4 patients",
# limits = c(3, 9),
expand = c(0, 0)) +
ggtitle("Fear of Bleeding: 20/100\nGREEN = MEDICATE") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
ggplotly(p1)p2 <- ggplot(tt,
aes(x = l2pb,
y = l2ps,
fill = lab)) +
geom_errorbar(aes(ymin = log2(10^4*(prob.stroke - 1.96*se.stroke)),
ymax = log2(10^4*(prob.stroke + 1.96*se.stroke))),
size = 0.5) +
geom_errorbar(aes(xmin = log2(10^4*(prob.bleed - 1.96*se.bleed)),
xmax = log2(10^4*(prob.bleed + 1.96*se.bleed))),
size = 0.5) +
geom_point(shape = 21,
aes(size = sqrt(N))) +
geom_abline(slope = m2$coefficients[2],
intercept = m2$coefficients[1],
linetype = "dashed") +
geom_ribbon(data = pi_20_80,
aes(x = l2pb,
ymin = upr,
ymax = 9.25),
fill = "green",
alpha = 0.2,
inherit.aes = FALSE) +
scale_x_continuous("Log2 of Bleeding per 10^4 patients",
expand = c(0, 0)) +
scale_y_continuous("Log2 of Stroke per 10^4 patients",
expand = c(0, 0)) +
ggtitle("Fear of Bleeding: 80/100\nGREEN = MEDICATE") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
ggplotly(p2)new_l2pb <- seq(from = 5,
to = 9,
length.out = 10*nrow(tt))
pi_20_80 <- data.table(predict(m2,
newdata = list(l2pb = new_l2pb),
interval = "prediction",
level = 0.6))
pi_20_80$l2pb <- new_l2pb
head(pi_20_80)p1 <- ggplot(tt,
aes(x = pb,
y = ps,
fill = lab)) +
geom_errorbar(aes(ymin = 10^4*(prob.stroke - 1.96*se.stroke),
ymax = 10^4*(prob.stroke + 1.96*se.stroke)),
size = 0.5) +
geom_errorbar(aes(xmin = 10^4*(prob.bleed - 1.96*se.bleed),
xmax = 10^4*(prob.bleed + 1.96*se.bleed)),
size = 0.5) +
geom_point(shape = 21,
aes(size = sqrt(N))) +
geom_ribbon(data = pi_20_80,
aes(x = 2^l2pb,
ymin = 2^lwr,
ymax = 400),
fill = "green",
alpha = 0.2,
inherit.aes = FALSE) +
geom_line(data = pi_20_80,
aes(x = 2^l2pb,
y = 2^fit),
linetype = "dashed",
inherit.aes = FALSE) +
scale_x_continuous("Bleeding per 10^4 patients",
limits = c(50, 400),
expand = c(0, 0)) +
scale_y_continuous("Stroke per 10^4 patients",
limits = c(0, 400),
expand = c(0, 0)) +
ggtitle("Fear of Bleeding: 80/100\nGREEN = MEDICATE") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
print(p1)p2 <- ggplot(tt,
aes(x = pb,
y = ps,
fill = lab)) +
geom_errorbar(aes(ymin = 10^4*(prob.stroke - 1.96*se.stroke),
ymax = 10^4*(prob.stroke + 1.96*se.stroke)),
size = 0.5) +
geom_errorbar(aes(xmin = 10^4*(prob.bleed - 1.96*se.bleed),
xmax = 10^4*(prob.bleed + 1.96*se.bleed)),
size = 0.5) +
geom_point(shape = 21,
aes(size = sqrt(N))) +
geom_ribbon(data = pi_20_80,
aes(x = 2^l2pb,
ymin = 2^upr,
ymax = 400),
fill = "green",
alpha = 0.2,
inherit.aes = FALSE) +
geom_line(data = pi_20_80,
aes(x = 2^l2pb,
y = 2^fit),
linetype = "dashed",
inherit.aes = FALSE) +
scale_x_continuous("Bleeding per 10^4 patients",
limits = c(50, 400),
expand = c(0, 0)) +
scale_y_continuous("Stroke per 10^4 patients",
limits = c(0, 400),
expand = c(0, 0)) +
ggtitle("Fear of Bleeding: 80/100\nGREEN = MEDICATE") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
print(p2)NOTE1: Deming regression requires standard errors for each point estimate. I do not have these estimates on the log scale. Hence, for the demonstration purposes, I’m showing regression models on the original scale.
NOTE2: how to do Deming regression prediction interval?
m1 <- lm(prob.stroke ~ prob.bleed,
data = tt)
summary(m1)
Call:
lm(formula = prob.stroke ~ prob.bleed, data = tt)
Residuals:
Min 1Q Median 3Q Max
-0.0060634 -0.0024344 -0.0001482 0.0016504 0.0059096
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.002773 0.002342 -1.184 0.263725
prob.bleed 0.688676 0.124877 5.515 0.000256 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.003809 on 10 degrees of freedom
Multiple R-squared: 0.7526, Adjusted R-squared: 0.7278
F-statistic: 30.41 on 1 and 10 DF, p-value: 0.0002564
m2 <- deming(prob.stroke ~ prob.bleed,
ystd = se.stroke,
xstd = se.bleed,
data = tt)
m2
Call:
deming(formula = prob.stroke ~ prob.bleed, data = tt, xstd = se.bleed, ystd = se.stroke)
n= 12
Coef se(coef) lower 0.95 upper 0.95
Intercept -0.003480583 0.001783199 -0.006975588 1.442144e-05
Slope 0.726042116 0.218915292 0.296976028 1.155108e+00
Scale= 2.014618
p3 <- ggplot(tt,
aes(x = prob.bleed,
y = prob.stroke,
fill = lab)) +
geom_errorbar(aes(ymin = prob.stroke - 1.96*se.stroke,
ymax = prob.stroke + 1.96*se.stroke),
size = 0.5) +
geom_errorbar(aes(xmin = prob.bleed - 1.96*se.bleed,
xmax = prob.bleed + 1.96*se.bleed),
size = 0.5) +
geom_point(shape = 21,
aes(size = sqrt(N))) +
geom_abline(intercept = m1$coefficients[1],
slope = m1$coefficients[2]) +
geom_abline(intercept = m2$coefficients[1],
slope = m2$coefficients[2],
linetype = "dashed") +
ggtitle("Least Squares (solid line) vs.\n Deming (dashed line) Regression") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
print(p3)sessionInfo()R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17763)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] MethComp_1.30.0 deming_1.4 plotly_4.9.3 ggplot2_3.3.3
[5] data.table_1.13.6
loaded via a namespace (and not attached):
[1] pillar_1.4.7 compiler_4.0.3 tools_4.0.3 boot_1.3-25
[5] digest_0.6.27 evaluate_0.14 jsonlite_1.7.2 lifecycle_0.2.0
[9] tibble_3.0.5 gtable_0.3.0 nlme_3.1-151 viridisLite_0.3.0
[13] lattice_0.20-41 pkgconfig_2.0.3 rlang_0.4.10 cli_2.2.0
[17] rstudioapi_0.13 crosstalk_1.1.1 yaml_2.2.1 xfun_0.20
[21] coda_0.19-4 withr_2.4.0 dplyr_1.0.3 httr_1.4.2
[25] knitr_1.30 generics_0.1.0 vctrs_0.3.6 htmlwidgets_1.5.3
[29] grid_4.0.3 tidyselect_1.1.0 glue_1.4.2 R6_2.5.0
[33] fansi_0.4.2 rmarkdown_2.6 farver_2.0.3 purrr_0.3.4
[37] tidyr_1.1.2 magrittr_2.0.1 scales_1.1.1 ellipsis_0.3.1
[41] htmltools_0.5.0 assertthat_0.2.1 colorspace_2.0-0 labeling_0.4.2
[45] lazyeval_0.2.2 munsell_0.5.0 crayon_1.3.4